Vibrational sidebands and dissipative tunneling in molecular transistors 
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Transport through molecular devices with strong coupling to a single vibrational mode is consid- 
ered in the case where the vibration is damped by coupling to the environment. We focus on the 
weak tunneling limit, for which a rate equation approach is valid. The role of the environment can 
be characterized by a frictional damping term S{uj) and corresponding frequency shift. We consider 
a molecule that is attached to a substrate, leading to frequency-dependent frictional damping of 
the single oscillator mode of the molecule, and compare it to a reference model with frequency- 
independent damping featuring a constant quality factor Q. For large values of Q, the transport is 
governed by tunneling between displaced oscillator states giving rise to the well-known series of the 
Frank-Condon steps, while at small Q, there is a crossover to the classical regime with an energy 
gap given by the classical displacement energy. Using realistic values for the elastic properties of 
the substrate and the size of the molecule, we calculate I-V curves and find qualitative agreement 
between our theory and recent experiments on Ceo single-molecule devices. 

PACS numbers: 73.23.Hk, 73.63.-b, 85.65.+h 
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I. INTRODUCTION 

In the emerging field of single-molecule electronics 
there is a large interest in transport through mesoscopic 
systems with strong electron-phonon coupling. There has 
been a number of experiments in which transport through 
a single molecule has been reported One exam- 
ple is the series of experiments by Park et al£ where it 
was shown that the current through a single Cqq molecule 
was strongly coupled to a single vibrational mode. The 
single phonon mode was associated with the motion of 
the molecule in the confining potential created by the 
van-der-Waals interaction with the electrodes. Later, 
similar devices with more complicated molecules were 
investigated;^ and they also showed excitation spectra 
which may be associated with emission of vibrational 
quanta. 

Theoretically, there has been a large amount of work on 
the problem of tunneling through a single level with cou- 
pling to phonon modes. In many experimental realiza- 
tions the tunnel coupling to the leads is rather weak, and 
the transport is dominated by the well-known Coulomb 
blockade effect. In this regime, where the transport is 
sequential, the use of a rate equation approach is appro- 
priate rather than a coherent scattering approach. Mo- 
tivated by the above mentioned single-molecule experi- 
ments, the rate equation approach has been used in a 
number of recent papersi^ai Some of these studies dealt 
with the issue of non-equilibrium phonon states, and the 
possibility of having negative differential conductance in 
such molecular systems:-'- 8 Physically, it is an essential 
question how the excited vibrational levels are allowed to 
relax, either through coupling to the environment, for ex- 
ample the phonons or plasmons of the metal substrate, or 
by virtue of the tunneling electrons^ In the case where 
the relaxation of the vibrational mode is faster than the 
tunneling rate one can assume an equilibrium phonon 



distribution. 

The coupling between the vibrational mode and the 
environment depends strongly on which vibrational mode 
is considered. For intra-molecular vibrations the lifetime 
can be very long^i^ However, for the experiments of 
Ref.S it was suggested that the vibrational motion was 
associated with a center of mass motion, which is coupled 
to the environment more strongly as we discuss in this 
paper. A sketch of the physical setup that we consider is 
shown in Fig. ^ 

In this paper, we assume that the tunneling rate is 
much smaller than the rate of relaxation to other degrees 
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FIG. 1: Illustration of the system considered in this paper. 
The molecule is attached to substrates, e.g., by van-der-Waals 
interactions, and the movement in this potential is modelled 
by springs with spring constants hu,\ and /cm, 2- When an 
electron hops onto the molecule, the force created by image 
charges or local electric fields causes a shift of the equilib- 
rium position of the oscillator and consequently emissions of 
quanta. The weights of the different final states are given 
by the well-known Frank-Condon overlap factors. The main 
objective of this paper is to consider the influence of damp- 
ing of the molecular motion by emission of phonons into the 
substrates. 
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of freedom. Hence, the usual rate equation approach is 
applicable, and we can assume that the phonons relax 
between each tunneling event according to a thermal bo- 
son distribution. We study a model of one single molec- 
ular orbital with strong Coulomb repulsion coupled to 
a dissipative environment. The dissipation is caused by 
coupling to phonon modes of the electrodes as well as 
electromagnetic modes and it is represented by a bath of 
harmonic oscillators. The description is thus similar to 
the well-known theory of Coulomb blockade in an electro- 
magnetic environment >i2ii^iiSii& However, there is a dif- 
ference in how the coupling to the environment appears. 
In the electromagnetic environment case, the tunneling 
of an electron results in a sudden displacement of the po- 
sition of the charge, while here the tunneling results in a 
sudden appearance of a force on the oscillator. For this 
reason, we go through the derivations in some detail and 
derive a general formula for the I-V curves. This gen- 
eral result does not depend on the nature of the environ- 
ment but we then specialize to two cases. We consider 
a molecule attached to a substrate, using a continuum 
model for the substrate, and compare our result to a ref- 
erence model featuring frequency-independent damping 
and quality factor. The I-V curves, as a function of the 
elastic parameters of the substrate and the size of the 
molecule, feature quite different line shapes as compared 
to the assumption of constant friction. 

To get a simple estimate of the importance of the 
coupling to the substrate, consider a model where the 
molecule position x is coupled to a one-dimensional sub- 
strate through a spring with spring constant km- For 
small substrate displacements, force balance gives that 
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(1) 



where u(z) is the substrate displacement, pro is the ID 
mass density, and v s is the sound velocity. At a given 
frequency u>, the outgoing soundwaves are u(z,lu) — 



, where a is a constant. Finding a from (1), we 



can insert it into the equation of motion for x to obtain 
the quality factor Q at the resonance frequency u>q: 
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where mo is the molecule mass with luq = kM / tjiq , and 
to = pidVs/luq is the mass of a wavelength long piece 
of the substrate. With realistic parameters for a 6*60 
molecule on a gold substrate, as was used in the experi- 
ments of Ref. |2j, the quality factor is between 1 and 10, 
and therefore we expect the broadening to be substan- 
tial. This furthermore confirms the assumption that, for 
this type of molecular device, relaxation through the en- 
vironment is much faster than through tunneling. 

The paper is organized as follows. The model Hamilto- 
nian is defined in Section [HI and in Section ITTT1 we derive 
an expression for the current from rate equations. The 
function that describes the tunneling density of states 



is then solved in absence of the dissipative environment 
in Section I1VI and with coupling to the environment in 
Section |Vj We discuss different models for the dissipa- 
tive coupling in Section IVII where we also discuss the 
physical implications. Section IVIII contains examples of 
I-V curves, and finally, a summary as well as a com- 
parison with the experiments of Ref. 2j can be found in 
Section ETTl 



II. MODEL HAMILTONIAN 

We consider a model of one single spin-degenerate 
molecular level coupled to two leads (generalization to 
more molecular levels is straightforward). The single 
level is coupled to the vibrational mode of the molecule 
through the charge on the dot. The coupling between 
the oscillator and the environment is included as a linear 
coupling to a bath of harmonic oscillators in the spirit of 
the theory by Caldeira and Leggett>i£ The model Hamil- 
tonian then reads 

H = H LR +H D +H B +H DB +H b&th +H Bhai . h +H T , (3) 

with 
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(4a) 
(4b) 

(4c) 
(4d) 

(4e) 
(4f) 



where the c£ , c ka and the d\, d a are creation and 
annihilation operators for the leads and the dot, respec- 
tively, x is the oscillator degree of freedom, {xj} de- 
scribes the set of environment degrees of freedom and 
uij , u>j their respective masses and frequencies, £0 is 
the onsite energy, and U is the Coulomb interaction on 
the molecule. The coupling constants between electron- 
oscillator and oscillator-bath are A and {f3j}, respectively. 
The lead electron energies are given by 



£,kct — £ ka Pa 



(5) 



where p a is the chemical potential of lead a. Finally, the 
tunnel Hamiltonian is 

H t= *fc«( c LcA + 4c fe< T,J- ( 6 ) 

ka. a—L.R 
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The tunneling amplitudes could in principle also depend 
on the oscillator position. In the experimental realiza- 
tions in Refs. yyjjj. this is probably a small effect since 
the oscillator length is of order a few pm whereas the 
tunneling matrix element changes on the scale of nm. 
For simplicity, we do not take any such non-linear ef- 
fects into account, but we note that a position depen- 
dence of the tunneling amplitudes, for example of the 
form exp(— xo/^t), could easily be included in the present 
formalism. 

The force acting on the charged molecule, represented 
by the term Hjjb, is caused by electric fields originat- 
ing from either static impurity charges or image charges. 
Since this force on the molecule is counteracted by a 
force on these charges and, hence, on the environment, 
we should in principle also include the interaction be- 
tween the environment coordinates and the charge on 
the molecule. This would in fact lead to a qualitatively 
different behavior since Ohmic dissipation is cut off at 
frequencies smaller than the inverse size of the total sys- 
tem, i.e., the inverse range of the interaction between 
the charge on the molecule and the charges in the en- 
vironment. This interesting subtlety was pointed out in 
Ref. [l8|. However, since the van-der-Waals interaction 
between molecule and the substrate is short range com- 
pared to the electrostatic forces, we will consider the force 
acting on the molecule as an external quantity, which is 
thus not coupled to the dissipative environment. We do 
note, however, that including such a coupling would in 
fact lead to a small discontinuity at the beginning of each 
step in the I-V curve. The experimental data of, e.g., 
Ref.H does not seem to suggest such a discontinuity, and 
we therefore specialize to the case where the environment 
coordinates are unaffected by the charge on the molecule. 

We now want to relate the coupling constants in the 
boson-bath coupling to the finite damping of the vibra- 
tional mode, which can be accomplished by studying the 
classical equations of motion. After removing the bath 
degrees of freedom, thereby neglecting the term H DB 
that will be removed by a unitray transformation below, 
we end up with the following equation of motion in the 
frequency domain: 



[uj 2 —Jq- 5(w)]x (w) = 0, 



where we have defined 
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(7) 



(8) 



which is complex in general and gives rise to frictional 
damping and a frequency shift of the bare frequency luq. 
In Section IVl Bl we will explicitly calculate S(u>) for the 
case of a molecule attached to a semi-infinite substrate. 

We eliminate the coupling term H DB of the Hamilto- 
nian by a unitary transformation similar to the one 
used in the independent boson model^ at the cost of in- 
troducing displacement operators in the tunneling term. 
However, since we are dealing with a somewhat more 



complicated system due to the coupling to the bosonic 
bath, the unitary transformation in Ref. has to be 
generalized. We define the transformation 

H = SHS\ S = e- iAn «, A = p l + y £ / p j l j , (9) 

3 

where n d = ^ a d\d a . Using that 



(10) 



it is a matter of simple algebra to show that the linear 
coupling term H DB cancels if we set 



1 = 
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to K 2 +5(0)]' > m^V 



(11) 



and the Hamiltonian then transforms into 

H = H LR + H D + H B + H hath + E Bhath + H T , (12) 
where 

Ht= £ t^ a (c{^J A d a + dU- iA c ka , a ) (13) 
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and 



H D = e J2 d Ua + Un dl n dV e Q = £ --\£. (14) 



Here, U = U — X£ is the Coulomb repulsion modified by 
the phonon mediated interaction. For a weak Coulomb 
interaction, this can result in a negative effective U, 
which was discussed in Ref. 123. 



III. RATE EQUATIONS AND CURRENT 
FORMULA 

We derive an expression for the current in the weak 
tunneling limit using the usual kinetic equation approach. 
As mentioned in the Introduction, the most important as- 
sumption here is that the tunneling rate is much smaller 
than all other time scales, which means that we can as- 
sume the vibrational degrees of freedom and the Fermi 
seas in the two electrodes to be in equilibrium at all times. 
For simplicity, we consider only two charge states and 
therefore let U = oo, which leaves us with only three 
states: empty, and occupied by either spin up or down. 
The probabilities for the three states are denoted Pq , P\ , 
and Pj , respectively. The rate equations are 



(15) 



which combined with the condition P + P, + P, = 1 has 
the solution 
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p = Pt = - 
1 T r 01 + 2i\ 



(16) 
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where Tio is the tunneling rate for tunneling from the 
empty state to a singly occupied state, and r i is the 
rate for the reverse process. Since the electron can tun- 
nel out of both left and right leads, both rates have left 
and right contributions: = Tfj + Tf-. The tunneling 
rates are calculated using Fermi's Golden Rule, thereby 
treating Ht of Eq. as the perturbation and assuming 
a thermal equilibrium distribution of the lead electrons 
and the phonon bath. Following standard derivations, 
we obtain 



1 10 — 1 o 



r a — r 
1 oi — 1 o 



dus 

duj 
2n' 



-F(u)n a (s +cj), (17a) 

>(-w)(l-n a (eo+w)), (17b) 
where we have defined the function 

F{uo)= dte tut F(t), F(t) = (e iA ^e- iA ), (18) 



in addition to the Fermi distributions of the two 
leads, n Q (e)=(e /3 ' e_eV °' 1 + 1) _1 , and the bare rates 
r a =27r^ fe |£fc Q | 2 <5(£fc). The function F has the proper- 
ties 

/DO J 
^-F{u>) = 1. (19) 

We can use Eq. (|T$)l to show that the rates in Eq. l(T7|) 

can be written as 



I\o — r a n Q , 
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(20) 



where we defined 



2^ 



F(u>)n a (uj + e ). 



(21) 



The current through the molecule is now given by 
7 = -e(2P rf -(P T +P i )r : 1 ) = 2 



pi? rL pi?. rL 
1 10 1 01 1 01 1 10 



r 01 + 2r 10 

Using Eq. (|20(l . this can also be written as 

2eT L T R fi R fi L (e^( E o- e ^) - e ^ e °- eV ^) 



(22) 



J = 



r L n L (2 + e P( e °- eV ^) + T R n R (2 + e^ e °- eV ^) 



(23) 



IV. WITHOUT COUPLING TO THE 
ENVIRONMENT 

We start by discussing the limit when the oscillator 
is not coupled to the environment, which means that 
thermal smearing dominates over dissipative broadening. 
However, we still assume that the coupling is stronger 
than the tunnel coupling so that the molecule equilibrates 
between each tunneling event. This section is thus equiv- 
alent to the results in other rate equation calculations, 





-2-10 1 2 



FIG. 2: Upper panels: current-voltage characteristics for a 
device without coupling to the environment for symmetric 
(Y R = r L ) and asymmetric (Y R — 0.05F L ) tunneling con- 
tacts. Lower panels: contour plot of the differential conduc- 
tance in the voltage-gate voltage plane, where eV g = eo- The 
curves have been calculated using the analytic result Eq. 12:j> . 
valid in the limit where the lifetime broadening of the os- 
cillator is negligible. The temperature is kT = 0.1cj for 
the thick lines and in the contour plots, while the thin lines 
are for kT = 0.025(^0- The bias is applied symmetrically 
Vl = V/2 — —Vr and in the I-V curves we take eo to be 0, 
1.5, and 2.5 times Uq. The current is measured in units of 
I N = eY L Y R (Y L +Y R ). 



but for completeness and later comparison we write down 
this limiting case. 

The phonon average is performed assuming thermal 
equilibrium, and we have 

F (t) = (e* , ^e-* , t°)'), (24) 
= exp {g (e-^* — l)(l + N)+g (e^* - l) N} , 

where 



1 ( I 



9= 2{T 



1 



m LO 



-, N = n B (u> ) ) (25) 



Here, g is an important parameter determined by the ra- 
tio of the classical displacement length and the quantum 
mechanical oscillator length. The evaluation of Fq(ui) 
from Eq. (|24|) is equivalent to the independent boson 
modelfS and using the result from there we get 

oo 

F (w)=2tt Pn(g)S(u-ncj ), (26) 

71 — — OG 

where 

P n (g) = exp(- 3 coth(&)) e ™ b /„ (— |^)) , b = 

(27) 

and is the modified Bessel function of the first kind. 
The finite-temperature result involves both positive and 



5 



negative values of n, corresponding to emission and ab- 
sorbtion of phonons, respectively. At zero tempera- 
ture, this reduces to having only positive values of n 
because of the factor e nt3uJ °^ 2 , and hence only emission 
processes are possible. In the limit T — > 0, we thus 
have a series of emission peaks at u> = ntoo for posi- 
tive n and with weights given by the Poisson distribution 
P n -> e- g g n /nl. 

The current can now be found from Eq. i|23[l . In Fig.|2J 
we show examples of current- voltage characteristics using 
Eq. ill'-il) for symmetric and asymmetric junctions. In the 
following, we study how the physics gets modified by the 
coupling to the environment. 



Solving this linear set of equations for the Green's func- 
tions and inserting the results into Eq. I|32[l . we obtain 



2guj a 



LO 2 — UJq — S(u) 
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(35) 



where we have defined S(u>) = S(uj) — S(0) and the 
experimentally observable renormalized frequency ljq = 
u) 2 +5(0). Using the function B(uj) thus follows as 



B(u) = -ig 



1 + n B {io) 



Im 



■S(w)] 



(36) 



V. WITH COUPLING TO THE ENVIRONMENT 

In presence of coupling to the environment, the eval- 
uation of the function F(t) in Eq. (|18|l is in principle 
straightforward since the Hamiltonian is quadratic in the 
oscillator and bath degrees of freedom. We obtain 

F(t)=exp(B(t)-B(0)), B(t) = (A(t)A(0)) , (28) 

where the operator A is defined in Eqs. © and iJUJ. The 
expectation value (. . . ) is to be evaluated with respect 
to H without the tunneling term. At this point, it is 
convenient to use the fluctuation-dissipation theorem, 



B(w) = -2lm[B R (uj)}(l +n B (u)), 



(29) 



to express B(t) in terms of the corresponding retarded 
Green's function 



B R (t) = -i9(t)({A(t),A(0)}) . 



(30) 



Here, ns(oj) = {er u — is the usual Bose function. In 
order to find this retarded correlation function, we define 
the following auxiliary Green's functions: 



G§{t) = -i8(t)j([O(jt),A(0)]), 

from which we obatin B R as 

B R = {tG R +Y J ^G R )l 

j 

P. 



(31) 
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(32) 



The equations of motion for these functions are in fre- 
quency domain given by 
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where now g — £ 2 /2£ 2 is defined with respect to the renor- 
malized frequency luq, i.e. £q — 1/too(I>o. This result can 
then be used to find 



F(t) = exp 



dio , 
2^ (e " 



l)B{w] 



(37) 



Eq. l|57|l is equivalent to the result for the Coulomb 
blockade of a single tunnel junction with coupling to the 
electromagnetic environment^*^ In the Coulomb block- 
ade problem, the tunneling density of states was related 
to the impedance as seen from the junction, here the 
I-V characteristic is in a similar way related to the fric- 
tional damping of the oscillator mode. In both cases, 
the low energy form of the spectrum is a power law at 
low temperatures. At small frequencies u) <S u) and zero 
temperature, we get the following power law behavior: 

F(a,)« u-\ a^hmf^MV (38) 

Furthermore, we use the trick by MinnhagenSi to find 
the F-function as the solution of the integral equation 

F(w) = - f ^F(0B(( - W )(C - w), (39) 

which is useful for the numerical evaluation of F . 



VI. MODELS FOR S(cu) 
A. Frequency-independent quality factor Q 

As a first attempt, we can assume S to be of the form 



S{lu) = S(0)+i 



LUUJq 



(40) 



(34) 



which leads to a frequency-independent quality factor Q 
of the single vibrational mode, similar to the Ohmic dis- 
sipation model of Caldeira and Leggett in Ref. ^3 The 
model is also similar to the Coulomb blockade problem 
of an LCR circuit^ which is described by the same 
formula. The limit Q — > oo is seen to coincide with 



6 



F(e)/2n 



(a) 


r- 1 






10 


J 








7.5 












5 












2.5 








1 J . . 








1 2 3 4 5 6 7 





1 2 3 4 5 6 7 




1 2 




(f) 













A 









2 



FIG. 3: The function F(£) and its integral for different values 
of g and frequency-independent Q. The curves have been 
calculated from Eq. IKffljl at zero temperature. We take Q = 
20, 10, 5, 2.5, 2 5 /tt, 0.1, 0.01, where g = 0.5 in (a) and (b), g = 
1 in (c) and (d), and g — 2 in (e) and (f). The curves have 
been displaced for clarity by multiples of (1,0) in (a), (c), and 
(e), and by multiples of (0.25,1) in (b), (d), and (f) (largest 
Q to the left). For large Q, the integrated function goes to 
the dissipationless results, where the step heights are given 
by the Poisson distribution (leftmost staircase in (a), (c), and 
(e)), whereas for small Q it goes to a step function at e = guj 
(vertical line). Note also that the differential conductance at 
the first step remains sharp while the higher order steps are 
smeared when Q > Q c . 



the results in Section IIVI since in this limit B(uS) — > 
2tt(1 + n B {uj)){±)5{uj ± u>q), and when this is inserted 
into Eq. (|3T|) . we get i|2*I)l. We also see that, for a crit- 
ical value of Q c = 2g/n, the function F and, hence, the 
differential conductance change from having a divergence 
at small energies to vanish at small energies. 

In Fig. |21 we plot the function F and its integral for 
different values of g and Q. It is clearly seen how the in- 
creasing dissipation smears the Frank-Condon steps. For 
strongly underdamped coupling to the environment, the 
steps are only weakly smeared, and still visible even for 
Q = 2.5. For the special value of 



Q = Q C 



2.9 



(41) 



the first step disappears and eventually for very small 
Q the function F goes towards a delta function, F — > 
27T 8 (u) — goJo). Physically, this means that in the small 
Q limit, the system relaxes immediately to the classical 
state and tunneling is only possible by paying the total 
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FIG. 4: The integral of the function F for g — 4 and different 
Q at zero temperature. We take Q — oo (dashed line), 20, 
8/tt, 0.1. At Q < Q c the steps disappear, while for Q < 
Q q the curves approach the classical limit, which is a step 
function at e = 4oj - The inset shows the F- function itself for 
the same parameters. 



classical energy cost of the displacement. To see this, 
we rewrite gui in terms of the coupling constant A and 
get X£/2 = \ 2 /2itioujq, which is the classical energy for 
displacing the oscillator by increasing the occupation rid 
by one. 

The crossover to the classical regime occurs when the 
lifetime of the oscillator, u> /Q, is comparable to the 
Heisenberg uncertainty time associated with the classical 
energy of the displaced oscillator, i.e. when (reinserting 
h) 



Zl = 

lu ~ X£/2 



=> Qq = ~ ■ 

9 



(42) 



The disappearance of the steps, which happens at Q cl 
is therefore different from the crossover to the classical 
regime. This is shown in Fig.^Jwhere we plot the integral 
of F for g — 4 for different values of Q. For Q = 20, 
the steps are only slightly broadened, while for Q = Q c , 
the steps are almost fully broadened but the line still 
follows the quantum behavior. Only for smaller Q do we 
approach the classical result, which is a step function at 

.9^0- 



B. Coupling to a substrate 

We consider a molecule of mass mo attached to a 
substrate that extends over the semi-infinite half-space 
z > 0. The case of two substrates, as is shown in Fig. ^ 
is a straightforward generalization and will be discussed 
at the end of the calculation. The 3D Lagrangian density 
for the substrate is given by»22 



1 



{d t uf -{vf -2v 2 t )[Vu 



(43) 
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where vi and Vt are the longitudinal and transverse sound 
velocities, and p is the mass density. This Lagrangian 
leads to the following equation of motion: 



dtU-v?V(Vu) + v 2 V x V x u = 



(44) 



Having in mind a small molecule attached to the ori- 
gin, the assumption of cylindrical symmetry around the 
z-axis seems reasonable. We define u r and u z as the dis- 
placements in radial direction and parallel to the z-axis, 
respectively. 

We consider the case that the molecule only exerts a 
total force J- perpendicular to the substrate surface: 



■ M 



2nr f '(r)u° z (r)dr 



(45) 



where u° z (r) is the parallel displacement at the surface 
defined by z — 0, and f(r) is a normalized distribution 
function, i.e. J 2nrf(r)dr — 1. This imposes the follow- 
ing boundary conditions on the stress tensor 



=o = 0, 



T zz \ z=0 = -Ff(r), 



(46) 



where the components T zr and T zz can be written as 
functions of the displacements u r and u z (see, e.g., 
Refs. |3 HH) • The solution is then a straightforward 
generalization of the procedure for a point source with 
f(r) cx 8(f) /t outlined by Lamb in Ref. |2^. We obtain 
in frequency space 



knfk 



Jo(kr)dk, 



H Jo G(k,uj) 
where we have defined the following quantities: 



(47) 



VtX 



k 2 - 



V 



if k 2 > 



t,i 



V 



t,l 



-it k 



<£■ ifr < 

>i,i v i,i 

.1 C\_ , r,,._\/,.2 



(48) 



G(k,u) = Ap L viv t k 2 - (\ l + 2hl)(v; + k 2 )vf 

+\ L k 2 v 2 +X L k 4 , (49) 

where /i£ and Aj, are the Lame coefficients, which are 
related to the sound velocities as 



vi = 



'\l + 2pl 



v t = 



(50) 



and fk is the Fourier-Bessel transform of the force distri- 
bution f(r), 



fk 



f(r)Jg(kr)rdr. 



(51) 



The (— )-sign in the definition of Vt,l m Eq. (|48|l is nec- 
essary for selecting the retarded response u) — > lu + irj 
corresponding to outgoing waves since the square-root 
function has a branch cut on the negative real axis. 



The total force T involves u° z (r) and vice versa, see 
Eqs. so that we obtain: 

r2nrf(r)u° z (r)dr = x (52) 
o 1 + R (u) 



where 



R(to) = k M 



v t Jo G ( k , w ) 



(53) 



The function of interest, S(w), can then be deduced from 
the equation of motion for the molecule in frequency 
space: 



Mw 2 io(w) = — km 



x (u>) 



2ttt f (r)u z (r)dr 



(54) 



Identifying ^m/wm with the bare frequency u>q, we ob- 
tain 



UJ — OJn + UJn 



R(u) 



x (lu) = 0, 



n TT~RM_ 

which implies upon comparison with Eq. (JJJ 



5(«) 



R{u) 



1 + R{u)' 



(55) 



(56) 



For the situation where the molecule is attached to two 
substrates it is simple to see that the function S(ui) be- 
comes instead 



S(u) = - 



kM,x 



^i(w) 



*CM,2 



R 2 (uj) 



M l + R x (u) M 1 + R 2 (u)' 



(57) 



where Ri,2 is given by Eq. I|53|l but with kja replaced 
by kMi,2 and, if the two substrates are different, with 
substrate parameters changed accordingly. However, be- 
cause of the lack of detailed knowledge about the actual 
geometry of the device, and since the coupling to the 
two sides of the junction is very likely to be asymmetric, 
we will in the following make the simplifying assumption 
that the molecule only couples to one substrate. 

Our results for R(u>) in Eq. (|53|l imply that the imag- 
inary part of S(tu) (which will eventually be responsible 
for the frictional damping) has contributions not only 
from extended waves in the substrate but also from waves 
that are confined to the surface, the so-called Rayleigh 
waves. Mathematically, this contribution arises from 
G(k, ui) being zero for a specific value of k. This value 
falls into the regime where both v t and v\ are real, i.e., 
where wavevectors k are larger than allowed for transver- 
sal and longitudinal waves, see Eq. fH^l^ 

In order to compare our result l|56|l to experimental 
data, we need to choose a specific model for the force dis- 
tribution function f(r). The most realistic model would 
involve a distribution in accord with the van-der-Waals 
potential, however, as a result of that fk is a rather in- 
volved function of k. For simplicity, we therefore choose 



f(r) = 



1 



2-kD 2 



c-^ D ,t.e. 



fk — 



1 



27iVl + k 2 D 2 ^ 



(58) 
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FIG. 5: Real part (dotted line) and imaginary part (solid 
line) of S(lu) as a function of frequency for a C&o molecule on 
a gold substrate and the particular choice D/Do = 1. The 
imaginary part tends towards zero linearly which is illustrated 
in the inset: the "quality factor" wc^o/ImS^) tends towards 
a constant at zero frequency. 



The parameter D is on the order of the width Dq of the 
molecule, e.g., Dq = 10. 4A for a Cqq molecule. For this 
model, we can explicitly extract 5(0), since then 



R(0) = 



oj 2 Ma 4 



64 (a 2 



■ 2 D' 



(59) 



where a = vi/vt- Note that i?(0) is proportional to the 



squared bare frequency Wg 
up with 

5(0) 



5(0) so that we end 



l + Q$R(p)/u& 



(60) 



This result has the particular effect on the damping co- 
efficient Im5(aj)/w that it is independent of D/Dq at 
zero frequency. We show plots of the real and imagi- 
nary parts of S(ui) in Fig. El The real part, and thus 
the renormalization of the bare frequency as a function 
of energy, goes to zero rather quickly, whereas the imag- 
inary part remains nonzero over a large frequency range. 
The latter is important for the damping since the quan- 
tity wa)o/Im5((jj) takes the place of the quality factor. 
The fact that this quantity tends towards a constant at 
zero frequency illustrates that the imaginary part of S(uj) 
rises linearly with u> for small cj. 

We plot the results for F and its integral in Fig.El First 
we note that the shape of the staircases is markedly differ- 
ent from the constant Q-factor model: they are asymmet- 
ric and less steep on the rising side with a rather sharp 
transition to the next step. The asymmetry is even more 
obvious in the peaks of F itself. We also note that a 
larger spread of the coupling over the surface, i. e. larger 
D/Dq, makes the peaks in F and the steps in its integral 
sharper and less asymmetric. For large D/Dq, the stair- 
case tends towards the large constant Q limit of before 
since then u /In\S(uj) grows rapidly with lu and at the 
same time 5(0) — > 0. 
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FIG. 6: The function F(£) and its integral for different val- 
ues of g, making use of 1561 for S(oj). The curves have been 
calculated from Eq. 13911 at zero temperature, assuming a Ceo 
molecule attached to a gold substrate. We take D / Do = 0.5, 
0.75, 1, 2, where g = 0.5 in (a) and (b), g = 1 in (c) and 
(d), and g = 2 in (e) and (f). The curves have been displaced 
for clarity by multiples of (1,0) in (a), (c), and (e) (largest 
D/Do to the right), and by multiples of (1.5,1.5) in (b), (d), 
and (f) (largest D/Dq at the bottom). The staircases in (a), 
(c), and (e) feature sharper and less asymmetric steps for 
larger D/Do but are clearly visible in any case. The asym- 
metry is even more apparent in the plots of F in (b), (d), 
and (f). In contrast to the constant Q-factor results in Fig. El 
even the first step in the staircase gets smeared for smaller 
D/Do- However, we recover the large constant-Q limit for 
large D/Dq. 



VII. I-V CURVES 

In this Section, we show a number of I-V curves 
using the expression in Eq. I|23[l at zero temperature 
based on the F-function, both for the case of frequency- 
independent quality factor and for the substrate model 
(|56H for S{lj) discussed in the previous section. 

In Fig. we show current-voltage characteristics for 
constant Q = 5 and g = 0.5, 1, 2. For this value of 
Q, the Frank-Condon steps are still visible. If we took 
even smaller values of Q (not shown) such that the steps 
disappear, the characteristics are still strongly modified 
by the electron-vibron coupling in the sense that a gap 
develops in the I-V curve. Such an effect was recently 
claimed to be observed in a different type of deviceiSl 

We also show I-V curves corresponding to a Cqo 
molecule coupled to a gold substrate, using the substrate 
model of Sec. lVI Bl see Fig- EI We display the I-V curves 
for g = 0.5, 1, 2, both for symmetric and for asym- 
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FIG. 7: Current-voltage characteristics for g — 0.5, 1, 2, 
and frequency-independent Q = 5 at zero temperature. In 
each panel, we have taken £o = 0, 0.5o;o, and u , and the 
voltage is applied symmetrically across the device, so that 
Vl ~ V/2 — —Vr. The current is measured in units of 
In = er L T a (T L + F R ). The panels on the left show the 
case of symmetric tunneling contacts, whereas the panels on 
the right side correspond to asymmetric tunneling contacts 
with Y L /V R = 0.05. 
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FIG. 8: Current-voltage characteristics for g — 0.5, 1, and 2 
using Eq. 156H for S(uj) at zero temperature, calculated for a 
Ceo molecule on gold. In each panel, we have taken £o = 0, 
0.5lj , and lu , and the voltage is applied symmetrically across 
the device, so that Vl = V/2 = —Vr. The current is measured 
in units of In = eT L T R (T L +F R ). The panels on the left show 
the case of symmetric tunneling contacts, whereas the panels 
on the right side correspond to asymmetric tunneling contacts 

with r L /r R = 0.05. 



metric tunneling contacts, however, we restrict ourselves 
to D / Dq=1 since the general features are very similar 
for other choices of D/Dq. Upon comparison with the 
frequency-independent Q-factor model, we note that the 
I-V staircases are in general less steep and smoother but 
still clearly exhibit the expected Frank-Condon steps. 



VIII. SUMMARY AND DISCUSSION 

A. Summary 

We have included broadening of the phonon sidebands 
due to frictional coupling of the oscillator mode within 
a kinetic equation approach. Since we have worked in 
the limit where the tunneling time is much smaller than 
the lifetime of the oscillator, we have assumed that the 
oscillator and the environment are in thermal equilibrium 
and in this case an analytical result for the current is 
obtained. 

In the reference model featuring the frequency inde- 
pendent oscillator quality factor Q, we recover the usual 
Frank-Condon physics for large values of Q. The transi- 
tion between the two different charge states is then given 
by the usual overlap of two displaced oscillator wave- 
functions, the governing parameter being the ratio of the 
displacement length I and the oscillator length £ Q , or g = 



£ 2 /2£q. For moderate quality factors Q > Q c = 2g/ir, 
the steps are smeared but still visible. For even smaller 
values of the quality factor, the decay time of the os- 
cillations becomes shorter than the quantum mechanical 
uncertainty time, which happens when Q < l/g. In this 
strongly damped case the tunneling process crosses over 
to a regime with a gap given by the classical displacement 
energy. 

We were also able to calculate S(uj) for a molecule that 
is attached to a substrate and showed how the molecule 
loses energy to the substrate. The model features similar 
general results as the constant Q-factor model, however, 
it is different in that the steps in the I-V curves rise 
more smoothly but feature a rather sharp transition to 
the next step, which then again rises up smoothly. The 
underlying reason is the peak structure of F, which ex- 
hibits asymmetric peaks due to the frequency dependent 
damping coefficient. We also note the dependence on the 
spread of the coupling over the substrate surface, param- 
eterized by D/Do, where our results tend towards the 
large constant Q limit for large D/Dq. 



B. Comparison with experiments 

We have also tried to fit thepresent theoretical re- 
sults to the experiments in Ref. |3, for which the theory 
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FIG. 9: Example of a fit to the experimental curves of Ref. 
using the substrate model 1561 for a Ceo molecule on gold, 
with g = 2 and D/Dq — 0.75. The dots are experimental 
data points for a gate voltage of 6.8V and positive bias volt- 
age, and the solid line is the theoretical curve. The smearing 
of the first step is seen to be reproduced well, while at the 
same time showing a sharp rise for the second step. This 
kind of smearing could not be produced by thermal smearing, 
which would smear both steps equally. However, it is not pos- 
sible to make consistent fits for the entire I-V curve and for 
different gate voltages. This suggests that the molecule might 
be changing position and/or coupling with changing voltages. 

should be appropriate since the tunneling broadening is 
much smaller than temperature, oscillator quantum, and 
observed widths. For these experiments, in which Cqq 



molecules were attached to two leads, it therefore seems 
likely that the broadening is dominated by coupling to 
the environment. 

A rough qualitative agreement, except for the steep- 
ness on the rising side of the steps, can be achieved for 
the frequency-independent quality factor if we assume 
g « 1 — 2 and Q « 2 — 6. However, in order to obtain 
quantitative agreement, it is necessary to assume differ- 
ent values for g and Q for different values of the gate and 
source-drain voltages. 

Our model for S(uj) that corresponds to a molecule at- 
tached to a substrate, see Sec. IVIBI features qualitative 
agreement with experiment if we assume g and D / D to 
be on the order of unity. The asymmetry in the peak 
structure of F actually provides for a better quantita- 
tive fit to the experimental data than is possible for the 
constant Q-factor model. This is illustrated in Fig. 
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